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ABSTRACT 



We explain the millisecond variability detected by Rossi X-ray Timing 
Explorer (RXTE) in the X-ray emission from a number of low mass X-ray binary 
systems (Sco X-1, 4U1728-34, 4U1608-522, 4U1636-536, 4U0614+091, 4U1735- 
44, 4U1820-30, GX5-1 and etc) in terms of dynamics of the centrifugal barrier, a 
hot boundary region surrounding a neutron star (NS). We demonstrate that this 
region may experience the relaxation oscillations, and that the displacements 
of a gas element both in radial and vertical directions occur at the same main 
frequency, of order of the local Keplerian frequency. We show the importance 
of the effect of a splitting of the main frequency produced by the Coriolis 
force in a rotating disk for the interpretation of a spacing between the QPO 
peaks. We estimate a magnitude of the splitting effect and present a simple 
formula for the whole spectrum of the split frequencies. It is interesting that 
the first three lowest-order overtones (corresponding to the azimuthal numbers 
m = 0, - 1, and - 2) fall in the range of 200-1200 Hz and match the kHz-QPO 
frequencies observed by RXTE. Similar phenomena should also occur in Black 
Hole (BH) systems, but, since the QPO frequency is inversely proportional to 
the mass of a compact object, the frequency of the centrifugal-barrier oscillations 
in the BH systems should be a factor of 5-10 lower than that for the NS systems. 
The X-ray spectrum formed in this region is a result of upscattering of a soft 
radiation (from a disk and a NS surface) off relatively hot electrons in the 
boundary layer. The typical size of the emission region should be 1-3 km, which 
is consistent with the time-lag measurements. We also briefly discuss some 
alternative QPO models, including a possibility of acoustic oscillations in the 
boundary layer, the proper stellar rotation, and g^-mode disk oscillations. 
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1. Introduction 

The process of accretion onto a compact object in the neutron star (NS) and black 
hole (BH) binaries has many remarkable similarities. Some part of accreting matter with 
large angular momentum forms a disk, while another part (participating in a sub-Keplerian 
rotation) undergoes practically a free fall accretion until the centrifugal barrier (CB) 
becomes sufficient to halt the flow (see e.g. Chakrabarti & Titarchuk 1995; hereafter 
CT95). Thus, one may expect that two distinct zones, a disk and a barrier, can be formed 
in the vicinity of a compact object which are likely to be responsible for the generation of a 
resulting spectrum. In this paper we consider a possibility of the dynamical adjustment of 
a Keplerian disk to the surface of a NS. 

The disk structure begins deflecting from a Keplerian one at a certain point to adjust 
itself to the boundary conditions at a NS surface (or at the last stable orbit, at Rq = 3Rs, 
in the case of a BH). The transition from a Keplerian to a sub-Keplerian flow may proceed 
smoothly, though it is very likely that a perfect adjustment never occurs. In general, the 
transition should take place through the setting up of the CB (where a centrifugal force 
slightly exceeds the gravitational force) within the adjustment radius. Here (see §2.1 for 
details) we discuss the formation of kinks and shocks in the supersonic regime of accretion 
flow as a possible physical reason for a super-Keplerian rotation. We suggest that in a 
region with a supcr-Kcplerian rotation a matter may experience the relaxation oscillations 
in the vertical and radial directions. These oscillations are expected to be in a resonance 
with the local angular velocity in a disk, and the variation in an emitting area caused by 
the oscillations around a transition point produces the QPOs in the X-ray flux. 

The principal observational consequence of our scenario is a correlation between the 
X-ray flux and QPO frequency. The relatively soft disk photons in the CB region are 
scattered off the hot electrons thus forming the Comptonized X-ray spectrum (Sunyaev 



- 5 - 



& Titarchuk 1980; hereafter ST80). The electron temperature is regulated by the supply 
of soft photons from a disk, which depends on the ratio of the energy release (accretion 
rate) in the disk and the energy release in the CB region (see Eq. [24] for more details). 
For example, the electron temperature is higher for lower accretion rates (see e.g. CT95, 
Figures 3, 4), while for a high accretion rate (of order of the Eddington one) the CB region 
cools down very efficiently due to Comptonization. 

In this paper we discuss the CB model as a possible mechanism for the ~ 1 kHz QPOs 
discovered by Rossi X-ray Timing Explorer (RXTE) in a number of LMXBs (sec Strohmayer 
et al. 1996, and Van der Klis et al. 1996, hereafter S96, and VK96, respectively; Zhang 
et al. 1996). These observations reveal a wealth of previously unknown high-frequency 
X-ray variabilities which are believed to be due to the processes occurring on, or in the very 
vicinity of an accreting NS. 

In §2 we schematically describe our CB model. We demonstrate that the main 
model parameter 7 (which is proportional to a ratio of accretion rate to viscosity and 
is in fact Reynolds number for the accretion flow) determines the QPO frequency. As 
the Reynolds number (7-parameter) increases, the CB moves toward a NS surface (or 
toward the inner edge of an accretion disk), and the X-ray flux at higher energies declines 
because of the cooling of the emission region (CB). In addition, as the Reynolds number 
(7-parameter) increases, the QPO frequency reaches its limiting value (§2.1). We estimate 
the characteristic frequency of the oscillations of the CB region in §2.2. In §2.3, we show 
that due to rotation of a disk the effect of the Coriolis force should result in the splitting 
of the main QPO frequency, and in §2.4, we present our estimate of the corresponding 
oscillation amplitude. We apply our CB model to the interpretation of RXTE observations 
of high-frequency QPOs in §3. In §4, we discuss the alternative possibilities including 
resonant acoustic waves in the hot boundary layer around the NS, stellar rotation. 
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disk nonradial oscillations for the QPOs with a frequency of ~300 Hz and g-mode disk 
oscillations. Finally, we summarize our conclusions in §5. 



Let us consider a Keplerian accretion disk around a BH or a weakly magnetized NS. It 
is generally believed that there is a transition layer in the vicinity of a compact object where 
an accreting matter adjusts itself either to the surface of a rotating NS, or to the innermost 
boundary of an accretion disk. In the transition layer the motion is not Keplerian, and is 
governed by the mechanism of the angular momentum loss by an accreting matter. Thus, 
we may define the transition layer as a region confined between the surface of a NS (either 
the last stable orbit or the corotating magnetosphere) and the first Keplerian orbit. In a 
Keplerian disk a matter is circularly orbiting with the angular velocity 



where G is the gravitational constant, M is the mass of a compact object, and R is the 
radius of an orbit. The radial motion in a disk is provided by the friction and angular 
momentum exchange between the adjacent layers that result in the loss of initial angular 
momentum by an accreting matter. The radial transport of the angular momentum in a 
disk can be written therefore in terms of the torque of viscous forces between the adjacent 
layers (see e.g. Shakura & Sunyaev 1973) 



2. Centrifugal barrier oscillations 



2.1. The formulation of the problem 




(1) 




Here M is the accretion rate, and Wr<^ is the component of a viscous stress tensor that can 



be expressed as 



(3) 
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where iJ^ is a half-thickness of a disk, and rj is the turbulent viscosity. The only parameter 
entering this equation is the ratio 

= ^ (A) 

which is nothing else but Reynolds number for the accretion flow. The viscosity determines 
the redistribution of the momentum in the flow, with both the particles and photons (in 
the radiation dominated region) participate in momentum transport between the shearing 
layers. Narayan (1992) has argued that the momentum transport (viscous diffusion) by the 
particles in the flow depends on the velocity of the flow if the sources supplying particles are 
"frozen" into the flow. Specifically, Narayan (1992) has demonstrated that if the particle 
sources comove with the flow at the bulk velocity V, and the particles diffuse with the 
velocity c^, then the effective diffusion (viscosity) coefficient scales as (1 — V'^/cl). Thus, 
the viscosity of the supersonic flow gets effectively suppressed provided that the collective 
plasma effects (Tsytovich 1977) do not have enough time to develop, which seems to be 
the case in the limit oi V Cs where the distribution of particles in the flow is essentially 
determined by the source (see Narayan 1992, Equation [3.5]). 

The very important consequence of the Narayan's result is that the regime of high 
7 (Reynolds number) can occur even for the low accretion rates, M, in the presence of 
high velocities of the flow. However, the radiative viscosity can be important under the 
radiation-pressure dominated conditions and in the presence of high velocities of the flow, 
which seems to be the case for the transition layer. The coefficient of radiative viscosity is 
given by 

V = ^nipUphc/, (5) 

where mp is the proton mass, Uph is the photon number density, c is the speed of light, and 
/ is the photon mean-free path. The photon number density is proportional to the accretion 
rate, M. Thus, it is important that if the KHz QPO phenomenon is mostly determined by 
the ratio of accretion rate to the radiative viscosity of the flow, M/rj, then we would expect 
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it to occur in the sources with both high and low accretion rates. 

We can calculate the distribution of angular velocity of matter in the transition layer 
by solving equation (2) with the appropriate boundary conditions. We adopt that at the 
inner boundary the angular velocity of matter either matches the rotation velocity of a NS 
(or velocity of the corotating magnetosphere at the magnetopause) , 

LU = ljq at R = Ro, (6a) 

or vanishes (e.g. in the case of a non-rotating BH), 

00 ^0 at R = Ro. (66) 

The outer boundary condition for the layer must require that the orbital motion of matter 
matches smoothly a Keplerian motion at some radius Rout- The latter means that the 
angular velocity of matter and its radial derivative should be equal to those for a Keplerian 
disk at the same radius, i.e. 

uj ^ujk at R = Rout, (7) 



and 



duj dcjk 



at R — Rout- (8) 



dr dr 

The boundary conditions (6)-(8) allow us to solve equation (2). Thus, we can unambiguously 
calculate a radial profile of the angular velocity in a transition layer a;(R), and derive the 
value of the outer radius. Rout- It is important, that the only parameter that controls the 
adjustment of the angular velocity in a transition layer to the Keplerian angular velocity in 
a disk is 7— parameter. 

Let us focus on how to get a formal solution that would describe a smooth transition of 
the Keplerian flow into the sub-Keplerian rotation. The corresponding problem for equation 
(2) should be formulated as the boundary problem for the second-order differential equation 



-9- 



with the three boundary conditions: one at the inner boundary, and two at the adjustment 
point, where the flow begins to deviate from the Keplerian motion in order to adjust itself 
to the sub-Keplerian rotation. Thus, we have three boundary conditions for the second 
order differential equation which implies that the position of the adjustment point can be 
determined uniquely (see e.g. Korn & Korn 1961, Ch. 9.3). Such a simple formulation 
is quite instructive: it guarantees the uniqueness of the solution and automatically 
fixes the position of the outer boundary. This is important because the position of the 
outer boundary of the transition layer is a priory not known. Due to this very general 
mathematical reasoning the specific radial dependence of the coefficients in the differential 
equation (2) is not of principal importance. 

To illustrate how one can fix the position of an adjustment radius, we shall present 
the solution of equation (2) for the case where the 7-parameter is a constant. Let us 
now introduce the dimensionless variables: angular velocity 9 — uj/ujq, radius r = R/Ro 
(Ro = 3Rs, Rs = 2GM/c^ is the Schwarzschild radius), and mass m = M/M©. We normalize 
the angular frequency by the value ujq = 2t: x 363 rad ■ s~^. The main reason for this is 
the remarkable 363-Hz QPO discovered during the type-I X-ray bursts from 4U 1728-34 
(VK96), and which is believed to be associated with the rotation frequency of a NS. 

In terms of these variables the Keplerian angular velocity reads 



The solution of equation (2) satisfying the boundary conditions (5)- (7) is 




(10) 



where 



'-'out — CNS^out 

-7 -2 
r . — r . 



and 



out — 



^K(rout)- 



out out 
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The value of the dimensionless outer radius, rgut, should be determined from equation (7), 
which can be now rewritten as 

^^out = Di7r-1 + 2(Cns - DOr„-^f (H) 

Thus, equations (10) and (11) represent the solution comprising both a BH (Cns — 0) and 
a NS (Cns — 1) cases. 

Editor, Please put Fig.l here 

In Figure 1 we plot a family of solutions of equation (2) for a NS case. The solid line 
shows a Keplerian solution described by formula (9), the dashed line shows the solution 
with a sub-Keplerian rotation in a transition region, and the dash-dotted line shows the 
solution with a super-Keplerian rotation in the upper half of a transition region. The 
formal explanation for these types of solutions is the fact that the occurrence of kinks in 
an accretion flow enhances (with respect to the Keplerian value) the absolute value of a 
radial derivative of the angular velocity at the outer boundary. The perfect adjustment of 
the Keplerian rotation to the sub-Keplerian inner boundary condition (see equation [10]) is 
possible only in a rather unrealistic case where at the adjustment point the angular velocity 
and its radial derivatives coincide with the corresponding Keplerian values. In this case one 
can uniquely determine the position of an adjustment point. However, if this continuous 
outer boundary condition (see [6] and [7]) is not fulfilled, then the solution of equation (2) 
subject to the inner sub-Keplerian boundary condition should necessarily have a regime 
corresponding to the super-Keplerian rotation. 

It is important, that the kinks generally imply the presence of the discontinuities in a 
radial derivative of an angular velocity of the flow. Thus, when this derivative exceeds the 
corresponding Keplerian value, the adjustment of a flow to the inner boundary condition 
most likely occurs through the formation of a region with a locally super-Keplerian rotation. 
The physical reason for a super-Keplerian rotation in a transition region is a partial outward 
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transport of the angular momentum (e.g. facilitated by the formation of kinks in a flow) 
due to viscous stresses. As a result of this redistribution of the angular momentum, the 
matter, before entering a regime of an essentially sub-Keplerian rotation (in the innermost 
part of a transition region), gets locally involved in a super-Keplerian rotation. 

Let us discuss the reasoning for the occurrence of a super-Keplerian rotation in more 
detail. It is well-known from the standard accretion theory (see e.g. Shakura & Sunyaev 
1973, Lipunov 1992, Popham & Narayan 1992) that, depending on the value of the viscosity 
(see Popham & Narayan 1992 for the detailed discussion of the viscosity effects), the flow 
may be either subsonic or supersonic. When the flow reaches a supersonic velocity the 
density of matter decreases (see Landau & Lifshitz 1988, equation [83,5]), and a cooling 
regime changes from the free-free cooling onto the Compton cooling (see §2.3, and also 
Zel'dovich & Shakura 1969). As a result, the temperature of matter increases, the effects of 
radiative conduction overtake the effects of hydrodynamic viscosity, and a strong isothermal 
shock develops (Landau & Lifshitz 1988, §95). Note also, that the change in the curvature 
of an accretion disk caused by an asymmetric illumination of a disk by a NS (Maloncy, 
Begelman & Pringle 1996) favors the formation of kinks. The kinks are characterized by the 
discontinuities in the coordinate derivatives of the velocity components. In their discussion 
of kinks in the hydrodynamic flows Landau & Lifshitz (1988, §96) emphasize a difference 
between the formation of shocks and kinks: the shocks are necessarily formed due to the 
speciflc continuous boundary conditions, while the kinks are always a result of singularities 
in the boundary or initial conditions. Thus, in the presence of kinks the adjustment of 
a Keplerian flow to the inner boundary condition can occur through a super-Keplerian 
motion. 

The Keplerian rotation can be adjusted to the condition at the inner boundary (see 
Figure 1) only if the absolute value of a radial derivative of angular velocity at the radius of 
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kink formation exceeds the corresponding Keplerian value. When a super-Keplerian motion 
occurs a matter piles up in the vertical direction thus disturbing the hydrostatic equilibrium. 
The restoring force due to the vertical component of the gravitational force prevents matter 
from further accumulation in a vertical direction and supports relaxation oscillations. The 
radiation drag force, which is proportional to the vertical velocity component, determines 
the characteristic decay time of the vertical oscillations. 

In our calculations illustrated in Figure 1 the dimensionless thickness of a region with 
a super-Keplerian rotation is about 0.2 (= 2.4 km, for a canonical NS of 1.4 Mq mass 
and 10 km radius). The interpretation of the ~ 27 /is time lags between the softer and 
harder components of a spectrum observed in kHz QPO from 4U 1636-53 (Vaughan et al. 
1997) in terms of the ST80 Comptonization model gives a similar estimate for a size of 
the hard emission region. The development of a zone with a super-Keplerian rotation may 
significantly affect the dynamics of a transition region: the CB may be set up on the way 
of a matter accreting onto a compact object. Another effect we shall address in this paper 
is a possibility of global oscillations of a ring-like configuration (see Figure 2 for an artistic 
concept of our model) formed by a matter accumulated near the CB. We will discuss these 
oscillations in Section 2.2 in a very schematic and semi-quantitative way. The detailed 
modeling of the transition region and analysis of its oscillation properties in the presence 
of the CB will be presented in a separate publication. We will also consider the effects of 
rotation on the oscillations in the case of an axisymmetric equilibrium disk structure. 

Editor, Please put Fig. 2 here 
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2.2. QPO frequency of the centrifugal bcirrier region 

Figure 1 shows that in the transition region there is a layer where the centrifugal force 
exceeds the gravitational force. This means that in this layer matter tends to pile up before 
it adjusts itself to a new hydrostatically equilibrium structure in the vertical direction. It is 
very likely that the relaxation to a new equilibrium structure will be accompanied by the 
nonradial oscillations of a layer. The detailed spectrum of oscillation modes depends on the 
physical parameters of the unperturbed configuration and is a subject of a special analysis 
(see e.g. Appendix B for a general introduction). In this paper, we would like to present 
an order-of-magnitude estimate based on a rather simplified model. Let us consider the 
small displacement, h(R, t), of the gas element in a vertical direction at the radius R. We 
assume that in hydrostatic equilibrium, the vertical gradient of the pressure P (= Pg + Prad) 
where Pg = 2pkT/mp is the gas pressure, and Prad = iTTpHdL/mp is the radiation pressure) 
is balanced by the vertical component of the gravity force Fg p(GM/R^)H(j(H(}/R). 
Thus, the equilibrium condition for a small vertical displacement reads (neglecting 
general- relativistic corrections) 

pGURl/R" = (TTpHdL/mp + 2pkT/mp, (12) 

where L(r) is the local radiation flux per cm^. Note that in this equation we have neglected 
the vertical component of a centrifugal force which is very small compared to the vertical 
component of a gravitational force. 

A small vertical displacement of a gas layer from the equilibrium will 
result in the restoring force fgr ~ — mpGMh/R^ and the radiation drag force 
|fr| ^ (4/3)£[GMHd)/R^](h/c). Here £ = L/L^d is the luminosity in Eddington 
units, and the inequality holds for a nonblackbody radiation field (the equality holds for the 
blackbody radiation field only). Thus, neglecting the density and pressure perturbations, 
we can write the equation describing the vertical oscillations of a layer (the equation of 
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motion) 

mpii = fgr - |fr|, (13) 

where dot means a time derivative. Now, by making use of equations (12) and (13), we can 
rewrite the oscillation equation in the form 

ii + ^[GMHd/cR=^]h + [GM/R=^]h = 0. (15) 

For a harmonic perturbing force (oc e'^*) the power spectrum of the oscillations reads 

where = GM/R^ is the eigen-frequency of the oscillator, and 7d = ^QgHd/c is the 
damping rate of oscillations. We can now estimate the oscillator's Q-value: 

Q = Qohd = c/(Hd£l^o) Z 300. (17) 

The emission temperature is determined by the accretion rate (see §2.3), and therefore the 
Q-value which scales as 

Q = c/(Hd£no) = l^^j J, (18) 

should also depend (through the temperature) on the accretion rate. Namely, the Q-value 
should increase with increasing a centroid frequency of the QPO, since the 
temperature of an oscillating layer decreases as the accretion rate increases. 

Editor, Please put Fig. 3 here 

In Figure 3 we present the QPO frequency as a function of the 7-parameter. This 
Figure implies that the QPO frequency should scale with the 7— parameter approximately 
as 7~°-^. Perhaps, this dependence manifests itself in the recent observational data on 4U 
0614+091 (Ford et al. 1997), indicating that the main QPO frequency depends on the 
blackbody flux. . 
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2.3. Temperature of a CB emission region 

The total count rate from the source increases with an increase in an accretion rate. 
This is accompanied by an increase in the observed QPO frequency (if the viscosity is not 
affected by the radiation, see §2.1). At the same time, the X-ray spectrum becomes softer, 
because of an increase in a supply of soft photons from the disk illuminating the emission 
region of the CB. 

The CB region can be treated as a potential wall at which the accreting matter releases 
its gravitational energy. This occurs in an optically thin region where the column density is 
of order of a few grams, or where the Thomson optical thickness Tq ~ a few. The amount 
of energy released per second is a fraction of the Eddington luminosity since the CB region 
is located in the very vicinity of a central object (~ 3-6 Rs). The heating of a gas due to 
the gravitational energy release should be balanced by the photon emission. For the high 
gas temperatures, Comptonization is the main cooling channel, and the heating of electrons 
is due to their Coulomb collisions with protons. Under such physical conditions the energy 
balance can be written as (see e.g. Zel'dovich & Shakura 1969, hereafter ZS69, equation 
[1.3]) 

F/To~Co-(^(a)£(T)Te/f(Te). (19) 

Here r is the current Thomson optical depth in the emission region (e.g. in a slab), a is the 
energy spectral index for a power- law component of the Comptonization spectrum, £(r) is 
a distribution function for the radiative energy density, f(Te) = 1 -|- 2.5(kTe/meC^), Te is a 
plasma temperature in K, Co — 20.2 cm is a dimensional constant, and, finally, 

(p{a) = 0.75q;(1 + a/3) if (p{a) < 1, otherwise (p{a) = 1. The latter formula is obtained by 
using the relationship between the zero- and first-order moments (with respect to energy) 
of the Comptonized radiation field (Sunyaev & Titarchuk 1985, §7.3, equation [30]). The 
distribution of the radiative energy density in the emission region s{t) can be obtained from 
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the solution of the diffusion equation (cf. ZS69, equation [1.4]), 

I d^g _ (F/c) 
3d^ ~ To ' 

subject to the two appropriate boundary conditions. 



(20) 



The first boundary condition must imply that there is no scattered radiation from the outer 
side of the emission region, i.e. 

^-^£ = forr^O. (21) 
dr 2 

In our case this condition holds at the inner surface of a slab facing a central object. We 
note that, in reality, one can expect some additional soft flux from a NS resulting in an 
additional illumination of the inner surface of a slab. In the following analysis we shall 
neglect this effect. 

The second boundary condition requires that at the outer surface of a slab the incoming 

flux should be equal to the external flux H (not to be confused with unrelated quantity 
which is a half-thickness of the disk) , i.e. 

1 d£ H 

-— = - atT = To. 22 
3dr c 

The solution of equations (20)-(22) provides us with the distribution function for the energy 
density 

<r) = ^{2 + 3ro[r/ro - 0.5(r/ro)^F/(F + H)]}. (23) 



c 

Thus, from equations (19) and (23) we get 

(^(a)Tero/f(Te) ^ 0.75 • 10^:^7^ K. (24) 

When H < F the spectral index a varies very little since a is a function of Tero/f(Tc) (see 
e.g. Titarchuk & Lyubarskij 1995 for the case tq ^ 2). Thus, as long as the external flux 
(due to the photons from the disk) is much smaller than the internal energy release (per 
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cm^ per second) in the CB region, the spectral index is insensitive to the accretion rate in 
the disk. The values of parameters consistent with equation (24) are: Tq ^ 5, Tc ^ 2 x 10^ 
K, and a ^ 1, which are characteristic of a hard state for the galactic BH and NS systems. 

When the external flux H becomes comparable to the internal energy release, F, the 
cooling becomes more efficient not only due to Comptonization, but also due to the free-free 
coohng, and therefore the electron temperature unavoidably decreases, Te <^ 10^ • F/(F + H) 
K (see CT95 for the numerical calculations of spectral indices and temperature). 

Editor, Please put Fig.4 here 

In Figure 4 we present the plot of a temperature in the CB emission region as a 
function of the H/F ratio. The latter is proportional to the mass accretion rate in a disk. 
We can calculate the temperature using expression (24), then we get the following formula 

for the spectral index (Sunyaev & Titarchuk 1980) 

a = (2.25 + 7c)'/' - 1.5, (25) 

where 7c = 7r^meC^/[3kTe(ro + 2/3)^f(Tc)] (not to be confused with unrelated quantity 7 
which is a Reynolds number of the accreting flow in the disk) . 

In the end of this paragraph we would like to present a rather elegant and instructive 
justification for the existence of a so-called "hot" solution for the case of a large optical 
depth (to » 1) of the energy release region. We must recall that a "hot" solution was 
discussed by TuroUa et al. (1994) and Zane, TuroUa & Treves (1997, hereafter ZTT97) in 
the context of a formation and structure of a static atmosphere around a non-magnetized 
steadily and spherically accreting NS. The statement regarding the existence of a "hot" 
solution (Turrola et al. 1994) was quite surprising, and this solution has been even 
interpreted as an artifact resulted from a perfect-reflection inner boundary condition 
adopted by the authors. Here we shall demonstrate that this solution follows directly from 
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our equations (19), (23), and (25) and assumption that H = OatT = To»l (the absence 
of an incoming flux from a boundary). 

After substitution of (23) into equation (19) and in the hmit of tq » 1 we get 

f2(Te) ^ 2.5. (26) 

This means that kTe = 120 keV which perfectly matches the temperatures presented in 
ZTT97 (see Figure 2, there). In our derivation of expression (26) we have used the relation 
a = 7c/3 (see equation [25] for a <^ 1), and we have calculated the energy density, e(r) 
(see equation [23]), for t/tq > 1/2. It is very interesting that in the limit of Tq ^ 1 the 
derivation of equation (33) is practically independent of an optical depth of the energy 
release region tq. Thus, we come up with exactly the same conclusion as in ZTT97: for a 
large optical depth of the energy release region tq the radiation and plasma are almost in an 
equilibrium, so that the plasma remains very hot ("hot" solution) because of the inefficient 
Compton cooling. 



2.4. Effect of rotational splitting 

In this Section we would hke to draw a reader's attention to the effect of rotational 
splitting of oscillation frequencies (Titarchuk & Muslimov 1997) that was ignored in all 
previous theoretical studies of QPOs. 

It is well-known (see e.g. Unno et al. 1979 and references therein) that in a rotating 
star the nonradial oscillations are split: (Jkim = c^ki + mQCki, where (7ki is the oscillation 
frequency in the non-rotating star, (Jkim is the frequency in the rotating star (in the 
co-rotating frame), m is the azimuthal number, Vt is the stellar rotation frequency, and Cki 
is an integral that depends on the stellar structure and on the eigenfunction. Thus, in a 
rotating star, a nonradial oscillation drifts at a rate mQCki relative to a fixed longitude of 
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a star. The effect of rotational splitting is due to the Coriohs force and is similar to the 
Zeeman effect in a magnetic field. 

The oscillations of a rotating disk must split as well, and the rotational splitting of the 
disk oscillations is particularly important for those modes whose frequencies are comparable 
to the frequency of a disk rotation relative to a distant observer. Our estimate (warranted 
when the relative correction to the eigenfrequency is small, a^^^ ja^ <^ 1) of the rotational 
splitting for a disk (see Appendix A, equation [All]) results in the following formula (see 
also Titarchuk & Muslimov 1997) 

2 



1 



n, (27) 



1 + STT^k^ + m^. 

where Qk,m is the oscillation frequency measured by a distant observer, is the 
eigenfrequency, is the local angular frequency of rotation of the oscillating region of 
a disk (as measured by a distant observer), s 1 is some function depending on a disk 
vertical structure (see Appendix A), and m and k are the azimuthal and vertical mode 
numbers, respectively. For the lowest-order modes with m = 0, — 1, and —2, the oscillation 
frequencies are 

^k,o = ^0, (28) 
^^1,-1 = -{^ - ^o) + TtA-T^' (29) 

-L ~|~ STT ~|~ J- 

and 

|Qi,_2| = fio + 2(Q - no) - , , V (30) 

i + STT + 4 

respectively. 



2.5. Oscillation amplitudes 



The observed amplitudes of the flux oscillations can be explained in terms of a variable 
area of the surface of a region emitting the hard radiation. The latter is produced by 
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the upscattering of soft radiation from a disk. The amount of soft radiation from a disk 
intercepted by the emission region, Lg, is proportional to its surface area, 27rHdRcb- The 
luminosity of hard radiation, Lh = A(ro,Te)Ls oc A(ro, Te)27rHdRcb, where A(To,Te), the 
enhancement factor due to Comptonization, is determined by the optical depth, tq, of the 
CB-emission region and by the electron temperature, Te (see e.g. Titarchuk 1994). Note 
that the CB oscillations do not affect the enhancement factor, since both the optical depth 
(or column density) and the temperature remain constant during the oscillations, so that 
their product ToTg ^ const (see equation [24]). We also point out that the radial and 
vertical displacements of a ring-like CB configuration should occur at the same frequency. 
Assuming that the volume of the CB configuration, Vcb = 27rHdRAR, does not change 
during the oscillations, it is straightforward to obtain the following relation between the 
vertical and radial displacements, SR^ and 5(AR), respectively, 

5{R) = S{AR) ^ -^AR. (31) 

To obtain an order-of-magnitude estimate of a relative amplitude of the luminosity 
variation due to the oscillations at the main frequency flo we can use the formula for a 
differential of the area of an emitting surface and equation (31), Scb = 27rRHd, 

5ScB = 27r(R - AR)5Hd ^ 27rR5Hd, (32) 

Then, for a relative change in the luminosity, which is proportional to a change in the area 
of an emitting surface, we get the following estimate 

5L/L ^ 5Hd/Hd. (33) 

Thus, even small variations of a height of a cylindrical CB area (5Hd ~ O.lHd can 
produce oscillations in the hard X-ray fiux at the level of SL ~ O.IL, which is of order of the 
observed ones. 
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We must note that an exact relation between the components of a displacement vector 
(see Appendix A for the notations) ^r, and depends also on the azimuthal and vertical 
mode numbers. For example, consider an individual oscillation mode with displacement 
^^^^ in an incompressible fluid. Then V • ^^^^ = 0, and we get a linear partial differential 
equation that determines a relation between C,r, iip, and C,z- Thus, in general, this relation is 
very complicated and should depend on the azimuthal and vertical mode numbers, m and 
k, respectively. 

The mechanism for the QPO emission we discuss in this paper implies that the 
QPO should have larger amplitudes in the hard tail of the spectrum that is produced by 
Comptonization in the CB region. This is in good agreement with the observations of 
the kHz QPO from 4U1636-536 (Zhang et al. 1996, fig. 2) for which the rms amplitude 
increases almost monotonically with the energy up to at least 20 keV. 



3. Theory and observations of kHz QPOs 



3.1. Main observational results 



In this Section we shall summarize the main observational results on the kHz QPOs 
that need to be understood. 

The launch of RXTE opened up a new era in the study of QPOs. Recently, the kHz 
QPOs have been discovered in the persistent fluxes of 8 LMXBs: Sco X-1 ( [van der Klis et al. 
[T996a| ), 4U1 728-34 QStrohmayer et at. 1996|) , 4U 1608-52 ( |van Paradijs et aL 1996| , [Berger 
\et al. 199q ), 4U 1636-53 ( [Zhang et al. 199q ), 4U 0614+091 ( [1-brd et al. 1996[) , 4U 1735-44 
( [Wijnands et al. 199q ), 4U 1820-30 ([Smale et al. 1996|) , and GX5-1 (jvan der Klis et al. 



[1996 b[ ). In addition, episodic and nearly coherent oscillations have been discovered during 
several type-I X-ray bursts from KS 1731-260 with a frequency of 363 Hz ( ^trohmayer ei 
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a 



I. 1996), during one type-I burst from KS 1731-260 with a frequency of 524 Hz ([Morgan 



fc Smith 1996|) , during 3 bursts from the vicinity of GRO J1744-28 with a frequency of 589 



Hz (IStrohmayer et al. 1996b| ), and during 4 type-I bursts from 4U 1636-53 with a frequency 



of 581 Hz ([Zhang et al. 1997a|) . The presence of a pair of the kHz QPOs is characteristic of 



almost all these observations. The centroid frequency of these QPOs ranges from 400 (4U 
0614+091, [l^brd et al. WS^ to 1171 Hz (4U 1636-53, [van der Khs et al. I996c0 . According 



to van der Klis et al. (1997) the observations of Sco X-1 also show the presence of additional 
two peaks at about 40 and 90 Hz simultaneously with the kHz QPO. 



For Sco X-1 ( [van der Klis et al. 1996(^ ), the difference in the centroid frequencies of two 



QPO peaks changes with time, whereas for 4U 1728-34 and 4U 0614+091 this difference is 
constant with time and does not depend on the count rate ( fStrohmayer et al. 1996| and 



Ford et al. 1996|) . In particular, for 4U 1728-34 this difference is always 363 Hz, the same 



as the frequency of nearly coherent oscillations observed during several bursts. For 4U 
0614+091, the difference in centroid frequencies for the two QPO peaks coincides with the 
centroid frequency of the third QPO peak which was observed during a 1/2-hour period. In 
general, the centroid frequencies of such QPOs range from 400 (4U 0614+091, [Ford et al 
to 1171 Hz (4U 1636-53, [van der Klis et al. l996c[ ). 



There are many similarities between the QPOs observed in different sources, and the 
most important common features are the following: 



1. For some sources ( Pord et al. 199"7| ) there is a correlation between the QPO 



centroid frequency and the count rate. However, Zhang et al. 1998 found that the 
frequency-count rate correlation is more complex in the case of Aql X-1 source. This 
correlation persists over a short time scale, from minutes to hours, and it apparently 
breaks dawn on longer time scale in this source (which has also been observed in other 
sources, e.g. in 4U 1608-52 [Berger et al. 1996]). However, in 4U 0614+091 ( [bbrd ei 
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al. 1997) the correlation is rather good over a long time scale. 

2. Large Q- values, up to ~ 10^. [There is also an indication that the higher frequency 
QPOs become more coherent as the frequency and the total count rate increase (|van| 
der Klis et al. 1997| ).] 

3. The rms amplitudes of QPOs range from the low values (at the threshold of 
detectability) to a maximum of 12 % in the RXTE/PCA band (2-60 keV). [For every 
source, where the data are available, the rms amplitudes show strong dependence on 
energy. For example, for 4U 1636-53 the rms amplitude at 3 keV is only 4%, while 
at 20 keV it is as high as 16% ( |Zhang et al. 1996| ). Furthermore, the variable hard 



X-ray flux seems to anticorrelate with the soft X-ray count rate and with a value of 
the highest QPO frequency (|Ford et al. 1996]) .] 



4. In Sco X-1, there is a good correlation between the frequencies of the 6 Hz, 40 Hz, 80 
Hz and the kHz QPO. Separation of two peaks anticorrelates with the kHz frequency 



and accretion rate, (see |van der Klis et al. 1997|) . 



5. The highest observed QPO frequencies fall in a remarkably narrow range from 1066 
to 1171 Hz ( [Zhang et al. 1997b| ). 



3.2. Interpretation of observational results 

We suggest that a set of peaks seen in the kHz QPOs is essentially due to the effect of 
rotational splitting of the main oscillation frequency (see Section 2.4). This effect should be 
very important for the accretion disk since the characteristic frequency of the gravitational 
oscillations of the CB region is of order of the Keplerian frequency. 

To estimate the split frequencies, we adopt that Qq ~ fl, which is a good 
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approximation for a nearly Keplerian accretion disk. Then, for the lowest-order 
modes with m = 0, — 1 and k = 1, 2, 3, 4 and 5, the oscillation frequencies are: 
fik,o ~ f^o, ^1,-1 ~ 2^0/(7:^8 + 2), ^2,-1 ~ 2^0/(47:^8 + 2), ^3 _i fa 2fio/(97r2s + 2), 
^4-1 2Qo/(167r^s + 2), and f^s -1 ~ 2Qo/(257r^s + 2), respectively. For m — —2 and k — 1 
we have -2I ~ ^^o[l — 4/(7r^s + 5)]. 

If we take = Qo/27r = 1200 Hz and assume that s ^ 0.7, we get the following 
frequencies (arranged in ascending order): u^-i ^ 14, ^ 20, 1^3-1 ~ 40, 1^2 -1 ~ 80, 
z/i _i ~ 270, and t^i,_2 ~ 800 Hz. These frequencies match the observed QPO frequencies in 
LMXBs. It is important that for a given mode with k > 2 the ratio between frequencies is 
practically independent of the function s and is solely determined by the quantum numbers. 
As has been pointed out by Hugh Van Horn (1997) the asymmetry between modes with 
m < (prograde modes) and modes with m > (retrograde modes) in our interpretation 
of QPOs in terms of the rotational sphtting may be akin to the result obtained by Carrol 
& Hansen (1982) that for the nonradial stellar oscillations the slow rotation enhances the 
stability of retrograde modes and makes prograde modes less stable. Note also, that the 
contribution of a (k, m) —component to the power spectrum is determined by a smoothness 
of a function over the corresponding coordinates. It is well-known from the Fourier analysis 
that the spectral power of the /-component scales as where n is an order of the highest 
existing derivative of a perturbation with respect to the corresponding coordinate (here 
I —m for azimuthal components, and I —k for z-components) . Thus, it is natural to 
expect the presence of only lowest-order azimuthal components (|m| = 0, 1, 2) in the power 
spectrum, since it is very likely that the perturbations along the 99— coordinate are rather 
smooth. On the contrary, the vertical perturbations should be essentially discontinuous, 
and their Fourier spectra should contain a large number of components with k > 1 whose 
contribution to the power spectrum ~ k"^. It is also worth noting that the spectral 
power of a (k, m)— component is related with the spectral power of a disk at the frequency 
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corresponding to this particular component. It is demonstrated (Lyubarskij 1996, Kazanas, 
Hua & Titarchuk 1997) that the density of the power spectrum for a Shakura-Sunyaev disk 
behaves as u!~^ at the low frequencies. Thus, because of this dependence, we can expect 
more power in the QPOs at lower frequencies. 

Therefore, we suggest that the rotational splitting of the main oscillation 
frequency of the CB region in the accretion disk may be responsible for the kHz 
QPOs in LMXBs. The observed three frequencies, ranging from 300 to 1170 
Hz, can be naturally interpreted in the following way: the lowest frequency of 
200-300 Hz corresponds to a mode with m = -1 and k = 1; the higher frequency 
of 800 Hz corresponds to a mode with m = -2 and k = 1, and the highest 
frequency is the main frequency of the oscillations. 

The main observational consequences of a CB model are: 

1. Correlation between the total count rate and the QPO frequency (see Figure 3). 
The latter increases with an increase in the accretion rate or, equivalently, with an 
increase in the total count rate (cf. Section 3.1, item 1) provided that the viscosity 
is determined by the turbulent motions. In the regime where the radiative viscosity 
dominates the frequency-count rate dependence may be more complex, and the kHz 
QPOs may even disappear. Zhang et al. (1998) show that in the source Aql X-1 the 
total X-ray flux (the total accretion rate) does not seem to be the only factor that 
controls the frequency of the kHz QPOs. We argue, that this frequency- count rate 
dependence is rather determined by Reynolds number ('-/— parameter). In other words, 
it is the ratio of the accretion rate to the viscosity that determines the frequency- count 
rate behaviour for the kHz QPOs. 

2. Correlation between the Q- value and the total X-ray flux. When the illumination of 
the CB region (the accretion rate) increases, the temperature of the CB region should 
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decrease (see Figure 4 and equations [18], [24]). 

3. Observed spectrum is a result of Comptonization of soft radiation from a disk in 
the CB region (ST80). The kHz oscillations is a manifestation of the mechanical 
oscillations of the hot CB region in the radial and vertical directions. 

4. Partial oscillation amplitudes should be larger at higher energies [exactly this 
correlation has been observed by [Zhang et al. 1996| (see Section 3.1 for details)]. The 
QPO frequency should also correlate with the total count rate and anticorrelate with 
the hardness of a spectrum (see Fig. 4, and also [I'brd et al. 1996| ). 



5. The CB effect should be more pronounced for the higher Reynolds numbers (7— 
parameters) for which the transition layer is closer to a NS surface, and for which the 
angular velocity in a disk and therefore the characteristic frequency of gravitational 
oscillations of the CB region i7o/2vr are larger. 

6. Spacing between the QPO frequencies and the main oscillation frequency (the highest 
QPO frequency) would decrease with increasing in the accretion rate (or increase of 
the main oscillation frequency). In fact the parameter s (see Eq. AS) monotonically 
increases from to 1 when the CB radius (or the oscillation) increases. The rotational 
splitting modes with the frequencies, 6 Hz, 40 Hz, 80 Hz, have to correlate with the 
main frequency 1200 Hz (see Eq. A9). 

7. The QPO oscillations are likely to become highly incoherent when the oscillating 
region plunges into the magnetopause, where the plasma turbulence enhances the 
viscosity and tends to suppress the oscillations. In this case the QPO frequency 
should be limited by the Keplerian frequency corresponding to the Alfven radius, 
Ra ^ 1.4 X 10'^(B/4 ■ 10^ G)^°-'^Mi^V-i/5 cm, where B is the magnetic field strength 
at the NS surface, Mig = M/(10^^ g s"^) and m = M/Mq. Note that the Alfven 
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radius and hence the highest frequency only slightly depend on the accretion rate. If 
the masses and values of the magnetic field strength for NSs in bursters are within 
a very narrow range, then one can expect that the highest frequencies for all these 
sources will also fall into a rather narrow range. 

8. Morgan, Remillard & Greiner (1997) reported the detection of 67 Hz QPO from GRS 
1915+105 by RXTE. They found the energy dependence of the 67 Hz QPO in their 
May 5 1996 observation, where the rms amplitude increased with the photon energy 
from a level of 1.5% below 5 keV to 6% above 13 keV. This correlation indicates that 
the 67 Hz QPO is associated with the most energetic spectral components visible in 
the PCA (Proportional Counter Array on board RXTE). In terms of our model it can 
be interpreted as the oscillation of the inner edge of the accretion disk (CB region), or 
as the gi— mode disk oscillation in the region of radius of 4 (see for details Appendix 
B). Such a correlation should be more pronounced in the model accounting of the 
flux variation in the converging-inflow produced hard tail (Titarchuk, Mastichiadis 
& Kylafls 1997, Titarchuk & Zannias, 1998), because of substantial variation in the 
exposure of the converging inflow region (the characteristic radius of which is 2rs) to 
the soft seed photons from the disk. 

The similar energy dependence is observed for 300 Hz QPO detected by RXTE 
(Remillard et al. 1996) from another superluminal source, GROJ1655-40. And thus 
we can argue that the QPOs observed from BH system detected only in the high energy 
hands reveal the scale of the high energy emission region and hence the nature of the 
high energy formation in the soft states of BH systems in terms of the converging 
inflow Comptonization. 



4. Alternative models 
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4.1. Acoustic waves 

The acoustic waves standing in the hot area surrounding a NS is an alternative 
possibihty for the interpretation of the kHz QPOs. The standing shocks above a NS surface 
can serve as boundaries for such waves. The existing controversy about the stabihty of 
such shocks will be probably resolved by 2-D simulations ( |Chen et al. 19971) indicating 
that after an initial transition phase of rapid shock movements a system relaxes to a 
stable configuration comprising the quasi-stable, hot torus around a NS. These simulations 
show that after passing a shock, the accreting gas becomes involved in the complex, 
predominantly tangential vortex motions with a very small radial Mach number in the 
postshock region. The resulting picture is quite distinct from a simple, quasi- virial spherical 
accretion: it allows for the existence of a hot coronal postshock region around a NS 
perfectly suitable for a sound wave propagation. The acoustic oscillations of such a region 
can manifest themselves as the kHz QPOs. Note that the NS systems seem to be more 
favorable, since the standing acoustic waves may oscillate between two surfaces, the stellar 
surface and the outer boundary. Moreover, these surfaces can reflect, absorb and even emit 
sound waves. 

The effects of an interaction of the shock surfaces with the sound waves are extensively 
discussed by Landau & Lifshitz (1989). In particular, the strong shocks are subject to 
corrugation instability when a shock surface becomes rippled. In the range of angles between 
the wave vector and a normal to the shock surface the reflection coefficient for an incident 
sound wave may exceed unity or even diverge. This means that the sound waves can be not 
only reflected but can also be spontaneously emitted by a shock (into the postshock region), 
so that a sound wave can receive energy from a shock. The damping of sound waves by a 
shock can also occur under certain circumstances. For example, for a sound wave falling 
normally from a postshock side onto a flat shock surface, the reflection coefficient is, in 
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general, < 1. Note that almost 100% reflection takes place when a gas is isothermal near the 
strong shock, and this may well be relevant to the physical situation under consideration. 
Due to Comptonization of the ionized gas the outer layers of a postshock region adjoining 
the shock are always isothermal. It is well-known that a similar situation occurs in the outer 
atmospheric layers of X-ray bursters, and in many other situations where the electrons in 
the outer, optically thin layers of gas clouds are heated by Comptonization and kept at a 
constant temperature equal to the spectral temperature of radiation emerging from a cloud. 
The ratio of specific heats is equal to unity and the refiection coefficient for a sound wave is 
almost equal to unity as well. 

However, it is unlikely that all sound waves to be perfectly refiected by a boundary 
of a hot region, regardless of its physical nature. The accretion disk, for example, may 
contribute to the wave damping, while the asymmetry may, as is mentioned above, result 
in a spontaneous emission and, therefore, in an enhancement of a sound wave. It is 
important that the standing acoustic waves in a postshock region is essentially a resonance 
response of a dynamical system (the lobe) to the background noise, which is always 
present in a postshock region and contains all possible frequencies. The balance between a 
continuous excitation of resonant modes and dissipation of sound waves determines their 
finite amplitudes. 

The waves under discussion are usual sound waves propagating in a fully ionized gas. 
The protons and electrons are coupled together by strong electrostatic forces, and the force 
driving the oscillations is a gradient of thermal pressure in an electron-proton mixture. In the 
situation characteristic of the inner boundary of an accretion disk the proton temperature, 
Tp ~ few keV is of order of the electron temperature, Tg. When entering a postshock 
region the protons instantly get much higher temperature, of order of the virial one, Tv ~ 
few MeV. However, because of frequent collisions with much cooler electrons the protons 
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lose their energy very fast, over a timescale tp^cooi ~ 5.8Tg/^(Mp/me)-'^/^/(neln A) s (Spitzer 
1956). Given Te ~ 10 kcV, ne - lO^^cm-^ and In A a; 13, we get tp,cooi ~ 2.5 x 10"^ s. The 
sound speed is therefore determined mostly by protons Cg = ((?p/9p)^/^ a; 0.98 x lO^Tg^^ cm 
s~^. Since protons and electrons are coupled by electrostatic forces, the frequent Coulomb 
collisions of protons determine the damping of sound oscillations, and we can use the proton 
mean-free path, Ip, in the standard expression for the wave damping in a fluid. The decay 
coefficient is 7d = jEml/E, where Em is an average rate of acoustic energy dissipation, and 
E is a total acoustic wave energy. Following Landau & Lifshitz (1989), we can write for the 
decay coefficient 7d ~ (27rz/)^(/p/cs), where u is a. sound frequency. 

Since the damping of sound oscillations in a hot plasma is mostly due to Coulomb 
collisions, the value Ip can be estimated as /p = (H/ro)(aT/o'c), where H and tq are the 
characteristic length-scale and the Thomson optical thickness of the plasma, respectively. 
The ratio of Coulomb to Thomson cross-sections is ctc/ctt ~ 2.35 • lO^Tg^, where T5 = T/5 
keV. The characteristic length-scale of an atmosphere, H, can be expressed in terms of the 
observed sound wave frequency as 

H=(7r/2)(c«/27rz/) = 2.4-Ty'//kHz-10^ cm, (34) 

where z/^hz = (^/l kHz). Assuming that H=const, we get u oc T^/^. The proton 
mean free path can be estimated (by using the above relations between Ip and H) as 
Ip — 0.1 ■ Tl^i/j^H^To"^ cm. Thus, one can find that 7d ~ 4 • 10~^uuizT^lTo^ , and the 
characteristic decay time t* = l/7d = 25To^'kHzT5^ s. This means that at high temperatures 
the oscillations should decay very fast. Note also that the oscillation frequency may vary 
as a result of a changing of the boundary condition at the shock. For example, in the case 
of a free boundary the oscillation frequency is a factor of 2 higher than that for a fixed 
boundary (the estimated size of a shock is a factor of 2 larger than the characteristic length 
scale given by [34]). The radiation pressure and/or viscosity can also play an important role 
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in damping of sound waves, despite the fact that the radiation pressure is small compared 
to a gas pressure at luminosities well below the Eddington limit (see §2.2). 

In the low (hard) state of a NS system the temperature of a postshock region is higher, 
~ few X 10 keV, and the acoustic waves die off rapidly. They can, therefore, exist only in 
the high (soft) state where Te ~ few keV, and the damping is small. 

The plasma oscillations usually considered in the laboratory plasmas do not propagate 
under the conditions we discuss here for a very simple reason. The effects of coUisionless 
plasma like plasma waves. Landau damping of ion-acoustic waves, etc. do not apply to 
our situation since the plasma we are dealing with is essentially coUisional. Given a short 
proton mean-free path, ~ 0.1 cm, and high velocity, Wp ~ 7 x 10^, cm s~^, the proton 
collision frequency is i^coi ~ 7 x 10^ Hz that exceeds a frequency of a kHz wave. 

As far as the observed amplitude of flux oscillations is concerned, it can be interpreted 
similarly to a CB model (Section 2), in terms of a variable area of an emitting surface. 

4.2. Rotation of a NS and oscillations of a disk 

Let us briefly address the issue of whether a NS itself may be a source of e.g. 363 Hz 
nearly coherent oscillations (S96). In general, in an oscillating system, the amplitudes of 
any two modes Aic^'^^* and A2e''^2*, may be either added to or multiplied by each other. 
When the excitation mechanism for the modes is the same, the addition rule applies, and 
Fourier spectral analysis will reveal two independent peaks in the power density spectrum 
(PDS), at uji and UJ2- This is the case, for instance, for a set of overtones of the same 
physical origin differing from each other by the mode quantum numbers only. When two 
oscillatory motions are of different origin, then we should multiply their amplitudes. The 
"beating" motion belongs to this category, including the "beat-frequency" model for the 
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low-frequency QPOs in LMXBs (see Lewin, Van Paradijs, & Taam 1993 for a review). In 
this case both (tui + UJ2) and \ijJi — uj2\ peaks should be present in the PDS of a system 
unless ijJi = UJ2- 

In LMXBs a NS is thought to be weakly magnetized and the lack of pulsations even in 
quite bright LMXBs apparently supports this idea. It is quite reasonable to assume that the 
nearly coherent ~363 Hz oscillations are associated with the stellar rotation: a relatively 
high temporal stability of a frequency and its independence of the spectrum and flux are 
the main arguments in favor of such a possibility. The 363 Hz oscillations occurred only 
during type-I X-ray bursts from 4U 1728-34 (other 3 LMXBs did not show any type-I bursts 
during the observations). The explanation this QPO in terms of a NS rotation encounters 
a difficulty that the emitting region should be seen only during a very short time (during 
bursts), which is, in principle, possible e.g. in a model of inhomogeneous combustion on 
a NS surface (see Bildsten 1995). The accretion disk trapped g'-modes is another possible 
mechanism that may be responsible for the stable a;3-oscillations (see Appendix B). Similar 
to a NS rotation, these oscillations, when superimposed on the intrinsic oscillations of a 
hot region around a NS may produce multiple QPO peaks. For the type-I bursts with 
a relatively short recurrence time (high values of M) the surface layers of a NS, remain 
convectively stable, and a combustion front propagates over the stellar surface very slowly, 
at Vcomb ~ 300 cm s~^ (Bildsten 1995). This means that a combustion inhomogeneity (e.g. 
associated with a more active burning in one zone where it started flrst) can persist for a 
long time, at least ~ 10^ s, provided that there is enough fuel. Thus, on a timescale of ~ 15 
s (typical duration of a type-I burst) a combustion inhomogeneity can manifest itself as a 
rather stable pattern of energy release on the stellar surface, which seems to agree with the 
4U1728-34 data. The narrow PDS peak imphes a high Q-value, again in agreement with 
the observations. 
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The period evolution of nearly coherent 363 Hz oscillations in 4U1728-34 during a burst 
can be understood in terms of the hydrodynamics of a "hot spot" produced by an accreted 
gas on the stellar surface (in a boundary layer). For example, as it has been suggested 
by Fujimoto (1988), the spreading of the accreted matter over the surface of a weakly 
magnetized NS accreting from a disk should be accompanied by a differential rotation in 
the very surface layers. This (radial) differential rotation gives rise to the hydrodynamic 
instabilities which will redistribute (through the development of a turbulence) the angular 
momentum in the surface layers of a NS. For the physical conditions characteristic of the 
surface layers of an accreting rapidly rotating NS the baroclinic instability (see Pedlosky 
1979, Knobloch & Spruit 1982, Zahn 1983, and Fujimoto 1987 for review) is perhaps the 
most important one. This instability is related to the fact that in a differentially rotating 
star in hydrostatic equilibrium the surfaces of constant density in general do not coincide 
with the surfaces of constant pressure. The growth rate for this instability can be estimated 
as (see e.g. Fujimoto 1988) r^^ ^ O.GOqcos^o, where ilo is the local angular velocity at 
the point with a colatitude 9o (< 7r/2). This estimate shows that the instability may 
develop as shortly as during 1-10 stellar revolutions. The development of the hydrodynamic 
instabilities means the occurrence of the largescale plasma flows and mixing of a matter 
in the surface layers of a NS that may eventually destroy a coherency of e.g. 363 Hz 
oscillations. 

Two other mechanisms have been proposed in the literature to explain the kHz QPOs. 
Recently, Klein et al. (1996) have proposed that these QPOs result from a turbulence 
occurring in the settling mounds of the neutron star polar caps. Their numerical simulations 
produce QPOs in the kHz range with fractional rms amplitudes of order of 1 %. Several 
QPO frequencies appear to be present in their simulations, with the highest ones are above 
2 kHz. Miller, Lamb and Psaltis (1997) have proposed a model in which the QPO with 
higher frequency is the Keplerian frequency at the sonic point beyond which the radial 
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inflow velocity becomes supersonic. The position of the sonic point is determined by the 
radiation forces which remove the angular momentum from the accreting matter. In this 
model the accretion flow is also modulated by the radiation flux from a NS with its spin 
frequency. This accretion-flow modulation causes the Keplerian frequency to beat with a 
NS spin frequency to produce QPO with a lower frequency. In this model, there is an upper 
limit for the QPO frequencies, either the Keplerian frequency at the stellar surface, or the 
frequency of the marginally stable orbit. 

We must admit that the variety of observational phenomena, such as e.g., the twin 
QPO peaks with the constant and (in some cases) variable difference between them, as well 
as appearance of a separate peak at the frequency corresponding to this difference need to 
be understood. We believe, the models of CB oscillations we attempted to address in this 
paper, may provide a reasonable and consistent interpretation for the most experimental 
facts related to the phenomena of kHz QPOs, both for the NS and BH binary systems. 

It is quite possible that ~ 200 — 400 Hz PDS peaks manifest either the proper spin 
frequency of a NS, or gi-mode oscillations in an accretion disk. The gi-modes can be observed 
at frequencies of ~ 200 - 500 and ~ (60 - 300) x (5/m) Hz for the NS and BH systems, 
respectively. The disadvantage of (7-modes from the observational point of view is that they 
should be localized in a rather small region of a disk with AR ~ GM/c^. 

The crucial test of potential observability of (?-modes in a disk around a relativistic 
compact object may be provided by the millisecond timing of BH transients in the high 
(soft) state where the luminosity is dominated by the thermal radiation from the accretion 
disk. Such 5(-modes should have frequencies ~ 200 • (5/m) Hz. The QPOs with the 
frequencies ~ 67 and ~ 300 Hz were recently detected by RXTE from GRS 1915+105 
(Morgan et al. 1997) and from GROJ1655-40 (Remillard et al. 1996), respectively, during 
the very high states of these sources. These observations may suggest the occurrence 
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of g'-mode disk oscillations, at least in these systems (see, also discussion in §3.2). The 
S'— modes can be observationally distinguished from the CB oscillations due to the fact that 
the gi-mode eigenfunctions are symmetric with respect to some radius in the disk (Nowak 
& Wagoner 1992), and thus the effect of rotational splitting will be smeared out for these 
modes (see Eq. Al), which is certainly not the case for the CB oscillations. 



5. Conclusion 

We propose a model that can self-consistently explain most of the observational facts 
on kHz QPOs. The main physical ingredients of our model are: 

1. The possibility of a super-Keplerian rotation (centrifugal barrier, CB) in a boundary 
region between a Keplerian accretion disk and a compact object. The formation of 
kinks, weak discontinuities in a supersonic accretion flow, makes an adjustment of the 
disk rotation to the innermost boundary (NS either, magnetosphere surface or the 
last stable orbit around a BH) to occur through the dynamical episodes of a locally 
super-Keplerian rotation. 

2. This dynamical adjustment is dictated by the turbulence of flow, Reynolds number, 
i.e. the ratio of the mass accretion and viscosity dictates the regime of the QPO 
behaviour. 

3. The excitation of relaxation nonradial oscillations in the CB region. 

4. The effect of rotational sphtting of the main oscillation frequency (of order of the 
local Keplerian frequency) of the CB region. 

We have demonstrated that the CB effect can explain various correlations and 
anticorrelations found in the observations of kHz QPOs (see §3.2 for details). Also, we 
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have shown that the sphtting of the main kHz frequency due to a disk rotation produces a 
discrete spectrum of frequencies that perfectly match the observed QPO frequencies (see 
§§2.4 and 3.2 for details). In our model the QPO frequencies are inversely proportional to 
the mass of a compact object. We have therefore suggested that a similar phenomenon 
should be observed in the BH systems too, with the QPO frequencies should be a factor of 
5-10 lower than in the case of the NS systems. In addition, in our model a typical size of 
the emission region is estimated to be ~ 1 — 3 km, in a good agreement with the time-lag 
measurements. We have briefly discussed some alternative models that may be relevant to 
the physics of QPOs including the standing acoustic waves in the postshock region around 
a NS, the effect of frequency modulation due to the proper rotation of a NS, and gi-mode 
intrinsic oscillations of an accretion disk. 

The authors thank NASA for support under grants NAS-5-32484 and RXTE Guest 
Observing Program. I.L. thanks the Isaac Newton Trust (Cambridge) and PPARC for 
support. The authors acknowledge discussions with Hugh Van Horn, Martin Rees, Guy 
Miller, Tod Strohmayer, Jean Swank, Nick White, Will Zhang. Particularly, we are grateful 
the anonymous referee whose constructive comments and criticism significantly improved 
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APPENDIX A. ROTATIONAL SPLITTING EFFECT 

In this Section we estimate the effect of a disk rotation on the QPO frequency. The 
detailed theory for a slowly rotating star is given by Ledoux & Walraven (1958) and Unno 
et al. (1979). Here we calculate the splitting of the eigenfrequency due to Coriolis force in 
the case of a disk geometry. 

It must be pointed out that in our analysis the correction to the main frequency 
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a^^^ produced by Coriolis force is always smaller than the main frequency ctq (even if 
the main oscillation frequency itself is comparable to the rotation frequency). Thus, to 
derive this correction, we are still allowed to exploit the linear approximation and use the 
unperturbed values for the displacement vectors. [The fundamental difficulty in the problem 
of nonradial oscillations in the presence of rotation is that the hydrodynamic equations are 
not t-invariant, as the Coriolis term is replaced by the corresponding term with the minus 
sign by time-inversion. This means that the eigenvalue problem differs from its complex 
conjugate, or, that the displacements ^ and ^* are not eigenvectors of the same eigenvalue. 
That is why the effect of Coriolis force on the nonradial oscillations can be taken into 
account only for the case of a slow rotation for which ^ and ^* can be approximated by the 
corresponding vectors determined from the eigenvalue problem without rotation, and 
^^^^*), respectively.] 

For the axisymmetric case (cf. Unno et al. 1979, equation [18.27]) the correction to the 
eigenfrequency reads 

^ " r AO) t(0)* ip> ' 

JD^b Sk,mSk,m*^^cb 

where ^^'^^(t, r, z) is the (k, m)-displacementent component. For a disk- (ring-) like 
geometry the displacement component should be calculated by using the (k, m) -harmonics 
of a complete set of eigenfunctions for a disk- (ring-) like configuration, {uk,m}- In equation 
(Al) the integration has to be taken over the disk area Deb of radius R, width AR, 
and half-thickness Ha. The appropriate displacement component ^^°m can be written in 
cyhndrical coordinates as (see e.g. Unno et al. 1979, equation [18.28]) 



Ukm, {A2) 



where the functions Uk,m satisfy the free-boundary conditions at z = Hd and z — — Hd (see 
e.g. Morse & Feshbach 1953), and the symmetry (with respect to the disk plane) condition. 
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These conditions are 

^(t, ^, Hd) = ^(t, ^, -Hd) = 0, {A3) 

and 

Uk,m(t, (f, Z) = Uk,m(t, (fi, -Z), (^4) 

respectively, where < z < Hd- We must note that the symmetry condition is not 
warranted in the case of a disk tilt above the plane z = on one side and below the plane 
on the other (Van Horn 1997). 

In the case of uniform rotation, the temporal and azimuthal dependences of 
eigenfunctions Uk,m are taken as e'^"^"^"^*^ (Unno et al. 1979). Thus, the eigenfunction 
component Uk,m can be written as 

Uk,m = e''"<^cos(7rkz/Hd)e-'^*, (^5) 

where m = ±1, ±2, k — 1, 2, 

The integral Ii in the numerator of equation (Al) is given by 

Ii = -27rHdRAR[21]mei^], (A6) 

where and are the average displacements over r- and </?— coordinates, respectively, for 
a disk- (ring-) like configuration. The integral I2 in the dominator of equation (Al) is 

I2 = 27rHdRAR [4^ + (^R/Ha) V^k^ + m^g] , (^7) 

where ^ is the average displacement over z-coordinate. Here 4, ^^p, and ^ are the 
average displacements for a disk configuration in radial, azimuthal, and vertical directions, 
respectively. To illustrate the effect of rotational splitting we may justifiably assume that 
ir ^ iifi-, then we arrive at 

(1) 2mQ 

^"l + STr^k^+m^' ^^^^ 
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where s(R/Hd) = (^zR/^Hd)^ is a function determined by the vertical structure of a disk. 
For more or less realistic structure of a disk s ^ 1 and depends on the ratio R/Hd, which is 
either ^ 1 or ^ 1. 

Thus, the frequency of the oscillations as seen by a distant observer is given by (cf. 
Unno et al. 1979| , equation [18.33]) 



1 - 



n. (A9) 



1 + STT^k^ + m^. 

In conclusion, we must note that this formula can also be obtained as follows. If 
the relative correction to the main frequency produced by the Coriolis term is small, 
then in the equation of motion for the perturbation we can use the displacement vector 
corresponding to the eigenvalue problem without rotation For a given oscillation 
mode the right-hand side of equation of motion (containing the gradient of a perturbed 
gravitational potential and the gradient of a perturbed gas pressure divided by density) 
can be replaced by the —al^^'^\ Then, by multiplying both sides of the equation of 
motion by ^^'^^* and integrating over the volume, we get the following quadratic equation 
cr'^ — 2a' a^^^ — (Tq = 0, where a' is the oscillation frequency in the corotating frame, and 
cr'-^'' is the correction to the main frequency due to the Coriolis term. The solution of this 
equation is a' = o-(^) + cro[l + (o-*-^Vo"o)T''^ ~ o-q + a'-^^ (if a^-^^ao < 1), which leads to the 
formula identical to (A9). Here we have ignored a solution with the minus sign, because it 
gives very high frequencies. 



APPENDIX B. ^-MODE DISK OSCILLATIONS 



Here we calculate the characteristic oscillation frequencies for the trapped g-modes 
in an accretion disk exploiting the results obtained by Nowak & Wagoner 1991, 1992, 
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and 1993 (hereafter NWl, NW2, and NW3, respectively). The two types of oscillations 
that are intrinsic to the accretion disks around the compact relativistic objects and are of 
potential interest for our present study are p- and gi-modes. The p-modes are longitudinal 
oscillations propagating along the disk and having displacements in the radial direction (in 
a disk plane). It is very unlikely that these oscillations can be observable. The g'-modes 
are mostly transverse oscillations, and the gravity is the main restoring force for them. 
Their wave vector is always in the plane of a disk, but the displacements are mostly in the 
perpendicular, z-direction. 

The essential feature of both p- and gf-modes is that they are trapped in a rather narrow 
range of radii. The possibility of the very existence of wave motions is essentially due to 
the characteristic feature of the Keplerian law, the stability of orbits. In the Newtonian 
gravity the orbits are stable at any orbital distance. The static general-relativistic effects 
e.g. in the case of Schwarzschild metrics result in that the stable orbits can exist only at 
R > 6GM/c^. For the modified Newtonian potential. 



the so-called epicyclic frequency (see NWl) — (c^/GM)(r~^ — 36r~^)^/^ vanishes at 
rm — 6, thus reflecting the lack of stable orbits at rm < 6. Here rm is the radial cyhndrical 
coordinate in units of GM/c^, and M is the mass of a central object. 



near this maximum vanishing near rm — 6 and at rm > rmax- The characteristic width of a 
trapping region is therefore Arm ~ 1, and the g^-modes die off exponentially outside this 
region over characteristic length scale ~ a disk thickness, h. The trapping at large radii is 
due to the pressure acting against the gravitational restoring force. The elasticity coefficient 
in the restoring force is ~ 1/R^, while the pressure is a much weaker function of R. As a 
result, the effective restoring force vanishes at rm — rmax 





The epicyclic frequency uj^ has a maximum at r, 



•\/60. The g'-modes are trapped 
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The equation governing the radial displacement Su as a function of r (NW2, equation 
[2.10]), 

u'clidn),. = -{uj' - ^Tn'){uj'' - <,J5u, {B2) 

coupled to the corresponding equation for the oscillations in z-direction (through a slowly 
varying with r function T[r]), fully describes the main propagation properties of waves, 
including a characteristic radial extent of the trapping region, and effective damping of 
waves out of this region. 

The dispersion relation for the waves describing both p- and g^-modes is 

u;'cy,^{u;'-^Tn'){u;'-u;l), (S3) 

where Q is a Keplerian frequency, and Cg is the sound speed. The values kr, T and fl in the 
dispersion relation are the values for the corresponding functions at rm ~ r^ax — 8 (NW2); 
'^e(rmax) — 0.029c^/(GM) , and T practically does not depend on r and varies from ~ 1 to 
~ 12 for the first five eigenmodes, depending on a mode number. 

The p-modes correspond to o;^ > tTQ^, u"^ > u^. As it has been pointed out in NW3, 
the p-modes are trapped in a very narrow region. Arm ^ O.lGM/c^, and, since they have 
displacements in the disk plane, these modes are unlikely to be observed. On the contrary, 
the gi- modes have frequencies below cUe and satisfy the conditions o;^ < cu^ and o;^ < ^TQ,^. 
Although the eigenvalue lj can be obtained from the exact solution of equation (B2) 
supplemented with a corresponding equation for the z-component alone, we may reliably 
estimate it from the dispersion relation (B3). For this purpose, let us rewrite (B3) in the 
form 

The frequencies o;^ < cv^ are allowed, and we can look for the solution of (B4) in the form 
uj = e ■ uje, where e^l. The hydrostatic equilibrium in the vertical direction implies that the 
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sound speed, Cg = 7?/^' = 7h^Q^ (SS73), satisfies the relation c^k^ = 7h^Q^(27r/A)^, where 
A is the radial wavelength approximately equal to the size (~ GM/c^) of a region with a 
mode trapping. Note that 7T[0(rinax)/t,;e(rmax)]Vs^ - (5/3)(l 12) • 2.25 • e'"^ > 1, and 
equation (B4) yields 

[l + T-i(27rh/A)2]~\ (S5) 

For a reasonable assumption that h ~ A(< 1), and given a typical value T ~ 6, we obtain 
e ~ 0.4. Thus, the value ~ a;e/2 is a typical value for a circular frequency for one of the 
lowest-order (yf-modes. Given fg = {uJc/2n){\ — 2/rjn)^^^ = 795/m Hz (here m = M/M©), 
these g'- modes should have frequencies f = uj/2t: ~ 200 — 500 Hz. These values perfectly 
fall in the range of frequencies between ~ 250 and ~ 360 Hz, observed from Sco X-1 and 
4U1728-34, respectively. 



Stability and damping of g-mode disk oscillations 

The Q- value for the ^f-mode oscillations is Q = f/Af, where Af is determined by 
damping. The damping of gi-modes is due to turbulent motions. The upward energy 
flux along the z-axis in g^-modes due to turbulent motions can be estimated as (Stein 
1967) Fgj,j^y (20 -=- 200)(4/Hd)^^Uo//t. A numerical factor in this formula depends 
on the assumed spectrum of turbulence. The energy density in a Keplerian disk is 
E — QY^^pjl — (^Vs/2)(R/Hd)^, and the damping constant can be estimated as 

In derivation of this formula we have used the same relation between the characteristic 
turbulent velocity, Uq, and the sound speed, Vg, as for the a-disks (SS73): Uq = aVg. 

Also, we assumed that the length scale of turbulent eddies 4 < Hd- Taking into account 
that Vs/Hd = v^/R = O and ^{x^ = 9 • 10^ • m s'^ , we get 7d = 4.3 • 10-3-0 s-\ 
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where ijj — m~^(Q;2)^[(4/-f^(i)i/3]^(7i)^/^(i?d/-R)^, and parameters 0:2, 4/Hd, and Ti 
are normalized by 10^^, 1/3, and 1 keV, respectively. The resulting estimate is 
Q~ 500 Hz/7d ~ 10^ • ip''^. The "experimental" Q-value for 4U1728-34 can be estimated 
from the fact that the oscillations persist for a whole duration of the burst, 15 s, which 
gives Qexp > 360 Hz/(1/15) = 5.4 • 10^. 



To estimate the g^-mode oscillation amplitude we can use the temperature distribution 
in a disk calculated by SS73. The relative contribution to the total luminosity, Ltot, 
from the annulus in a disk confined by the radii from —7.5 to 8.5, where the g^-modes 
oscillations are trapped, can be estimated as AL/Ltot — 0.04. For a disk isentropic 
in z-direction, as discussed in NWl-3, the formula PV = const (7 = 5/3) gives a 
relation between the temperature and thickness of a disk: Th'''"^ = const. Thus, 
|5T/T| = (2/3)5h/h. The flux from the radii between 7.5 and 8.5 will oscillate with the 
amplitude (5(AL)/AL = 45T/T = (8/3)(5Hd/Hd, and we can write 



which means taht ~ few % is an upper limit for the flux oscillations due to g^-modes in a 
disk. 



Luminosity variation produced by a disk gr-mode oscillations 
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Fig. 1. — Dimensionless angular velocity, 9, versus a dimensionless radius r in the disk. A 
family of solutions given by expression (10) (see §2.1) is presented for the case of a NS. The 
solid line shows the Keplerian solution described by formula (9). The dashed and dash-dotted 
lines correspond to the adjustment radius r=: 1.3 and 1.39, respectively. In our calculation 
the 7-parameter (see equation [8]) is equal to 15. 



Fig. 2. — Schematic diagram illustrating the centrifugal barrier (CB) model. The CB 
oscillates in vertical, radial and azimuthal directions around the vicinity of the adjustment 
radius (see the text for details) . Some fraction of the soft photons from a disk are radiated 
in the observer's direction, while most of the photons from a disk illuminate the CB region 
and get upscattered there to higher energies. 



Fig. 3. — Plot of f X m (where f is the QPO frequency in Hz, and m = M/M© is the 
dimensionless mass of a compact object) versus 7-parameter (oc M). 



Fig. 4. — Temperature of the CB region as a function of ratio H/F (where H is the external 
flux of soft photons from a disk that illuminate the CB region, and F is the intrinsic energy 
flux from the CB region. We performed our calculations for the Thomson optical depth of 
the CB region tq = 5. 
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